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We study the spreading of damage in the one-dimensional Ising model by means of the stochastic 
dynamics resulting from coupling the system and its replica by a family of algorithms that interpolate 
between the heat bath and the Hinrichsen-Domany algorithms. At high temperatures the dynamics 
is exactly mapped into de Domany-Kinzel probabilistic cellular automaton. Using a mean-field 
approximation and Monte Carlo simulations we find the critical line that separates the phase where 
the damage spreads and the one where it does not. 

PACS numbers: 



I. INTRODUCTION 

It is well known that most models studied in equilib- 
rium statistical mechanics, such as the Ising model, are 
defined in a static way through the equilibrium Gibbs 
probability distribution associated to the Hamiltonian 
of the model. It is desirable from the theoretical and 
numerical point of view to assign a dynamics to such 
models. The stochastic dynamics introduced by Glauber 
is the prototype example of a dynamics assigned to 
a static-defined model. The numerous versions of the 
Monte Carlo method [2j, used in statistical mechanics 
are also examples of dynamics assigned to static-defined 
models. All of them are markovian processes that have 
the Gibbs probability distribution as the stationary dis- 
tribution. In general they are either a continuous time 
process governed by a master equation ^ l a. [j , [a[ p,Ll B or 
a probabilistic cellular automaton 0, l9l llfl TnTll2L Il3| . 
The latter is defined by a stochastic matrix, whose el- 
ements are the transition probabilities, and the former 
by the evolution matrix, whose nondiagonal elements are 
the transition rates. 

If we wish to simulate, for instance, the Ising model 
we have to choose one of the possible stochastic dynam- 
ics since there are many. Having decided which dynamics 
to use, that is, having decided which probabilistic rules 
to use, we realize that there are several ways of doing the 
actual simulation corresponding to the chosen probabilis- 
tic rules. For instance, for the case of the probabilistic 
cellular automaton used by Derrida and Weisbuch ^(| 
to simulate the Ising model, and which will concern us 
here, there are several ways of realizing the dynamics. 
We may use the so called heat-bath algorithm [lj] or the 
algorithm introduced more recently by Hinrichsen and 
Domany |l5j or any other we may invent. These algo- 
rithms govern the movement of the system in phase space 
and they may be called stochastic equations of motion in 
phase space. Different algorithms may be the realization 
of the same probabilistic rule or stochastic dynamics. 

The description of a system either by the equation of 
motion or by the time evolution of the probability are 
equivalent. An analogy can be made with the Brownian 



motion which can be described either by the Langevin 
e quat ion or by its associated Fokker-Planck equation P, 
13. lid IT^ I . The first is a stochastic equation of motion of 
a representative point in phase space whereas the second 
governs the time evolution of the probability distribution 
in phase space. 

In the study of spreading of damage 0, 0, 0, 0, 
EHH E3, |H |H Hi it has been realized that algo- 
rithms that are realization of the same probalistic rules 
may yield different results for the spreading of damage 
EE El El El, and they usually do. The spreading 
of damage is a procedure through which we may study 
the sensibility of the time evolution of systems with re- 
spect to the initial conditions. The procedure amounts 
to couple the system with a replica of it, each of them 
following the same equation of motion. The coupling is 
acomplished by the use of the same sequence of random 
numbers. The equation of motion for each system to- 
gether with the use of the same random number define 
the equation of motion for the coupled system from which 
we obtain the joint transition probability [f3 . Il3| for the 
coupled system. 

Suppose one uses an algorithm to couple a system and 
its replica. This will lead us to certain joint transition 
probability. If another algorithm is used, which is also a 
realization of the same transition probability for a single 
system, the joint transition probability will be distinct. 
The correlation between system and replica will also be 
distinct and, in particular, the Hamming distance which 
is a measure of the damage spreading will be different. 
For example, in onc-dimcnsional Ising model, the heat- 
bath (HB) algorithm will give no spreading of dam- 
age whereas the Hinrichsen-Domany (HD) algorithm E3 
will exhibit a spreading of damage above a certain tem- 
perature This is an expressive example that the 
spreading of damage is not a intrinsic static property of 
a given system, but depends on the algorithm, or the 
stochastic equation of motion, we use to perform the ac- 
tual simulation [TH |22| . 

In this paper we introduce a family of algorithms, or 
equations of motion, spanned by a parameter that inter- 
polates between the HB and HD algorithms. The asso- 
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ciated transition probability corresponds, for all values 
of the parameter, to the Derrida-Weisbush (DW) proba- 
bilistic cellular automaton If we use this family of 
algorithms to study the spreading of damage, as we will 
do here, the parameter will have no effect on each system 
separately since for any possible value of the parameter 
the algorithm is related to the same transition probabil- 
ity. However, the joint transition probability will depend 
on the parameter and the properties of the system, in- 
cluding the damage spreading, will also depend on the 
parameter. 

A remarkable property of the dynamics introduced 
here is that at infinite temperature it is exactly mapped 
into the Domany-Kinzel (DK) probabilistic cellular au- 
tomaton El. This gives support to a conjecture by Grass- 
berger [22| according to which the generic class of damage 
spreading transitions is the same as the directed perco- 
lation to which belong the transition ocurring in the DK 
probabilistic cellular automaton. 

II. SINGLE SYSTEM 

Let us consider a one dimensional lattice where at each 
site one attaches an Ising variable Oi that takes the values 
+1 or —1 and denote by a = {<ii} the set of all variables 
of the lattice. The time evolution of the probability Pi(cr) 
of state a at discrete time i is given by 

Pl+l{<T')=Y J W{(T'\(T)Pt{cT) (1) 

a 

where IT(er'|<7) is the transition probability from state a 
to state a' which, for a probabilistic cellular automaton 
is given by [8| 

W(a'\a) =Y[w PC A((T' i \a) (2) 

i 

where wpcA(o'' i \a) is the probability that site i will be in 
state a[ in the next step given that the present state of the 
system is a. The DW probabilistic cellular automaton 
hOl for the one dimensional Ising model is defined by 

WpcA^'iW) = WDw{e'i\vi-i,Gi+i) (3) 

with 

iu.Dw(+l|0i-i,0i+i) =Pi{<r) (4) 

and 

u>rw(-l|0i-i,Oi + i) = 1 -Pi(a) (5) 

where 

g-/3J(<Ti_l+<T, + l) 

The site i assumes the state +1 with a probability Pi(a) 
that does not depend on the central site i. If we choose 



TABLE I: Transition probabilities for the DW probabilistic 
cclular automaton 
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the linear size of the system to be even the dynamics is 
decomposed into two independent dynamics for each sub- 
lattice. It is possible to show 01 that the DW probabilis- 
tic cellular automaton has as the stationary probability 
distribution the Gibbs probability distribuiton associated 
to the Ising model, namely, 

P(a) = |exp{/3jJ> j( x l+1 } (7) 

i 

where f3 = so that it defines a stochastic dynam- 

ics that can be assigned to the Ising model. 

The transition probabilities WDwiPiWi-ii are 
shown in Table I where we used the parameter p defined 
by 

P - e 2pJ + e -2/3J W 

The actual computer realization of a probabilistic cel- 
lular automaton can be made in several ways. Here, we 
introduce a family of algorithms that are possible real- 
izations of the DW probabilistic cellular automaton. It 
has a free parameter a that interpolates between the HD 
and HB algorithms. At each time step all sites of the lat- 
tice are updated in a synchronous way by means of the 
following algorithm, or equation of motion for the spin 
variables, 

<?i = sign{pi(o-) - (9) 
if Vi-x = a l+ i and 

a[ = sign{(a - &)(1 - a - - 6)} (10) 

if <7;_i 7^ Ci+i where & is a random number identically 
distributed in the interval [0, 1]. 

When a = one recovers the HD algorithm 0] 

er- = sign{pj(er) - £;} er;_i = er. 1+ i (11) 

o* = -signl^ ~ &} 0i-l 7^ Oi+i ( 12 ) 

and when a = 1/2 one recovers the HB algorithm 0, Hfij 

a[ = sign{pi(o-) - &} (13) 

It is straightforward to show that the algorithm defined 
by Eqs. © and (|10|) yields the one-site transition prob- 
ability given by Eqs. (0J and JSJ for any value of the 
parameter a. 
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III. COUPLED SYSTEM 

Let us denote by a — {c,:} and r — {t{\ the config- 
urations of the system and its replica, respectively. All 
sites of the system and its replica are updated in a syn- 
chronous way according to the algorithm 



if (Ji—i = <Ji+i and 



if ^ a l+1 and 



t[ = signer) 



if Tj_i = r i+1 and 



sign{(a-^)(l-a-^)(--^)} 



(14) 



(15) 



(16) 



(17) 



if T i-i 7^ T i+i- Notice that the random number is the 
same for both systems. 

The coupled system will be described by a four-state 
probabilistic cellular automaton defined by the time evo- 
lution 



of the joint probability Pi (er; r) of state (<r; r) at discrete 
time I where W(a'\ r'\a; r) is the joint transition proba- 
bility from state (a; r) to (ct';t'), and given by 

W((t';t'\(j;t) = Y[w(a' i ;Tl\a i -i,a i+1 ;T i -i,n +1 ) (19) 



^From the stochastic equation of motion given by Eqs. 
(fH|) . (fT~S|> . l(T^|) . and l(T7|) . we deduce the joint transition 
probabilities i«(<7^; r/lci-i, <7i+i; Tj_i, t^+i) that the site i 
of the system and the replica assume the values a[ and t[, 
respectively. The resultant joint transition probabilities 
are displayed in Table II and are valid for < a < p. 
For p < a < 1, the algorithm yields a joint transition 
probability which is independent of a and is the one that 
results by formally replacing, in Table II, a by p. The 
joint transition probability fulfill the following properties 

w(a' i ; T^\ai-i, a i+ i; n-i, r i+1 ) = WDw(&'iWi-i,Vi+i) 

(20) 

2J w(a' i ; t[ |(Tj_i, (T i+ i; Tj_i,T i+ i) = w D w(t~1 T i+ i) 

a' 

( 21 ) 

which contemplates the condition that the system and 
the replica follow their own dynamics independent of the 
coupling. 
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tem 



II: Joint transition probabilities for the coupled sys- 
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The joint transition probabilities satisfy also the fol- 
lowing properties, (a) Reflection symmetry in which the 
states of sites i — 1 and i + 1 are interchanged, that is, 
<-> (Ti+i and Ti-i <-> Tj+i. (b) System-replica sym- 
metry in which the states of the system and the replica 
are interchanged, that is, <Ji «-> Tj for all sites, (c) Up- 
down symmetry defined by the transformation <-> — cr.; 
and Ti <-> — Tj for all sites. 

The Hamming distance, that characterizes the spread- 
ing of damage, is defined by 



(22) 



which is also the order parameter related to the damage 
spreading phase transition. 



IV. RELATION WITH THE DK AUTOMATON 

In this section we show an exact relation between the 
stochastic dynamics defined in Section III and the DK 
probabilistic cellular automaton . If we let rji be the oc- 
cupation variable attached to site i, that is, r\i — or 1 ac- 
cording to whether site i is empty or occupied by one par- 
ticle, then the transition probabilities u^if (?7i|?7t-i, Vi+i) 
of the DK cellular automaton is given by 



w DK (1|0Q) = 



w DK {l\Ql) = w DK {l\lQ) = Pl 



W DK (l\ll) = p 2 



(23) 



(24) 



(25) 
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The DK cellular automaton displays a critical line in the 
phase diagram pi versus P2 that separates the absorbing 
state, for which the density of particles is zero, and the 
active state, for which the density is nonzero. 

Now, let us denote by rji the coupling variable associ- 
ated to the dynamics of Section III that takes the value 
1 or according whether ai ^ Ti or at — Ti respectively, 
given by 



1 



(1 - am) 



(26) 



The relation between the Hamming distance and the cou- 
pling variables is just 



* = (Vi) 



(27) 



The joint transition probabilities in the variables rji and 
<7j are defined by 



= w(a' i ;T' i \a i ^ i ,a i+ -i;Ti^ 1 ,T i+ i) 



(28) 



where r,- = cr,(l - 2jj,-) 

Summing over the coupling variable we get the follow- 
ing property 

w(a' i \r)' i \(j i -i,(Ji + i\'qi-i,rii + i) = WDw^'Mi-i^i+i) 

(29) 

which contemplates, as in Eq. H2U|) . the condition that 
the system follows its own dynamics independent of the 
coupling. The main property we wish to show, however, 
is that for infinite temperature, that is, for p = 1/2 we 
have 

with the DK transition probabilites defined by P2 = 
and pi = 1 — 2a. This means that the subsystem defined 
by the variables {rji} follows a dynamics identical to the 
DK probabilistic cellular automaton. ^From relation l|27|l 
it follows that the Hamming distance coincides with the 
order parameter of the active state displayed by the DK 
automaton. 

Yet for the case p = 1/2, it is easy to show that the 
joint transition probability satisfies the property 



= w DW (er - \<ji-i , <j i+ i)w DK {r)'i\r)i-i , r\ %+ \ ) 



(31) 



with the DK transition probabilites defined by P2 = 
and pi = 1 — 2a. Therefore, the c-subsystem and the 
77-subsystem are statistically independent. 



V. MEAN-FIELD SOLUTION 

Dynamic mean-field approximation has already been 
used to study systems in nonequilibirum stationary states 
[L3 u3 ■ Here we set up equations for an approxi- 
mate solution of the equation that governs the time evo- 
lution of the coupled system. We start by writing down 
the equations that give the time evolution of the one-site 
and two-site probabilities, namely, 

Pz + x(oi\ti) = ^2 X! w( <7 i; r ih I ^;To,T 2 ) 



(32) 



and 



f5?+l(01,<73;Tl,T3) = ^ ^ w(cri;Ti|cr ,er2;To,T 2 ) 



(70, 0-2.CT4 r 0: T 2,T4 



X w(<7 3 ;T 3 \a 2 ,<T4;T 2 ,T4:)P£(cro 1 (T2,a4 : ;T Q ,T2,T4) (33) 

^From now on we will drop the subscript £ and use the 
prime superscript for quantities calculated at time I + 1. 
To get a set of closed equations we use the approximation 

P{&Qi 02, 04; t 0i r 2 , t 4 ) = 



P(cr , 02; t , T 2 )P(a 2 , 04; t 2 , r 4 ) (34) 



P(0 2 ;r 2 ) 



which defines the dynamic mean-field pair approxima- 
tion. 

The probabilities P(ai; Ti) and P(<7i, 03; Ti, T3) cannot 
be considered all independent variables. Taking into ac- 
count that they should have the reflection symmetry and 
the system-replica symmetry and, in addition, assuming 
the up-down symmetry the probabilities are related as 
follows 



P(-; + ) = P(+;-) = -tf (35) 



P(+- 



P(-;-) = P(+;+) = -fi (36) 



P( ; ) = P(++;++) = A (37) 



-) = P(-+;— ) = P(— ;+-) 
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P(- 



D 



P(-+;-+) = P(+-;+-) = C 



(38) 



(39) 
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D 



E 



(40) 



(41) 



These seven variables are not yet independent. Only 
three of them can be considered independent which we 
choose to be >F, B, and D. The others are related to 
them by the relations 



n = 2P{- 



A = P{++) - 2B - D 



C= l --P{++)-U> + D 



(42) 



(43) 



(44) 




FIG. 1: Phase diagram in the plane a versus p where p is relate 
to temperature by ■ The continuous line corresponds to the 
mean-field approximantion and the circles to the Monte Carlo 
simulations. 



E= -2B-D 
2 



(45) 



where P(+) and P(++) are the one-site and two-site 
probabilities corresponding to a single system. The exact 
solution of the one-dimensional Ising model gives P(+) = 
1/2 and P{++) = [1 + (tanh/3J) 2 ]/4. 

^From the time evolution given by Eqs. (|32fl and i|33|) 
and using Eqs. (|43[1 . 144fl . and l|45|) we get the following 
closed equations for 'F, D, and B 



= 2 7 D + 8aB 



(46) 
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d2 d2 tj td 

B' = 2aB — 4a 2 — — 4a 2 2a 7 (48) 

* <F v ' 



where 



and 



7 = 1 - 2p 



a = — + p — 2a 



(49) 



(50) 



A stationary solution of the evolution equation is such 
that the stationary probability P(&] t) is zero when a ^ 
r which corresponds to no damage spreading (\1/ = 0). 
/From Eqs. i|4T)|) . l|T7|) . and (0HJ| we may obtain solutions 
with damage spreading ( , F ^ 0). The transition line is 
obtained by a linear analysis of stability of the solution 
around *F = and by assuming that the variables B and 
D vanishes linearly with 'F. Taking the limit 'F — > we 
find a transition line given by the implicit equation 

(1 - a) 2 7 3 - 4a(3a 2 - 5a + 2) 7 2 + 



+ 4a 2 (13a 2 - 16a + 5) 7 - 8a(3a - 2) = (51) 

whose solution is shown in the phase diagram of Fig. 
1. In particular, when a — (corresponding to the HD 
algorithm) we have 7 = 2(1 — a) which substituted in the 
equation for the transition line gives 



1 - 9a + 33a 2 - 59a J + 53a 4 - 20a 5 = 



(52) 



whose solution is a = 0.696173 from which we get p = 
0.196173 so that J/k B T c = 0.352597 and T c = 2.83610. 
When a =1/2 (correspoding to the HB algorithm) there 
is no transition. 

At infinite temperature, p = 1/2, the mean-field tran- 
sition line gives a = 1/6. Now, using the relation 
pi = 1 — 2a obtained from the equivalence with the DK 
automaton, and taking into account the result p\ = 2/3 
obtained in [12| in the pair approximation for the DK au- 
tomaton we have a = 1/6 in coincidence with our present 
result. 



VI. NUMERICAL SIMULATIONS 

Our numerical simulations resulted in the transition 
line shown in Fig. 1. When a = we have ob- 
tained p = 0.285(1) which gives J/k B T c = 0.230(1) and 
T c = 4.35(2) i n ag reement with the result by Hinrichscn 
and Domany [jjl, namely J/ksTc — 0.2305. At infi- 
nite temperature, p — 1/2, the numerical results give 
a transition at a = 0.0955(1). Now, using the relation 
Pi = 1 — a obtained from the mapping of our model into 
the DK cellular automaton, we obtain p\ c = 0.809(1) in 
agreement with previous Monte Carlo numerical results, 
namely p lc = 0.8095(5) 28]. 

The determination of the critical line was obtained by 
using the time dependent method 0, H2, . We started 
with two one-dimensional lattices (system and replica) 
with L = 1000 sites. Both lattices were initialized with 
completly indepedendent random configurations so that 
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FIG. 2: Time dependent Monte Carlo simulations for the 
damage survival probability P for a lattice with linear size 
L = 1000. Numerical data are shown for a — 0.075 and 
p = 0.450, 0.453, 0.455, 0.457, and 0.460 from bottom to top. 



half the spins were damaged at the begining ("3/ = 1/2). 
The update was done in a synchronized way by using 
the algorithm defined by Eqs. (HJ|, 03, (HD, and (fTTjl . 
with the same random number for both lattices. The 
damage surviving probabilites P(t), obtained by taking 
the averages over 2000 samples, were collected from t = 1 
to t = 1500 Monte Carlo steps. At the critical point we 
expect the following asymptotic time behavior 

P{t) ~ t~ s (53) 

Therefore, a double-log plot of P versus t will be linear 
at the critical point. In Fig. 2 we show how the critical 
value of p was found when a — 0.075. Several values of 
p, the ones shown in Fig. 2, were checked in order to find 
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